package ua.lviv.franko.solvers.oneDimention;

import java.util.ArrayList;
import java.util.Iterator;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;

import ua.lviv.franko.FEMShapeFunction1D.IFEM1D;
import ua.lviv.franko.integration.ConstantFunction1D;
import ua.lviv.franko.integration.FunctionWithoutDerivatives1D;
import ua.lviv.franko.integration.GaussQuadrature;
import ua.lviv.franko.integration.IFunction1D;
import ua.lviv.franko.integration.Quadrature;
import ua.lviv.franko.matrix.RowType;
import ua.lviv.franko.matrix.hash.lu.HashMatrixLUStyle;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.mlsFunctions.oneDimension.BaseFunctionCoupled;
import ua.lviv.franko.mlsFunctions.oneDimension.ICalculate;
import ua.lviv.franko.mlsFunctions.oneDimension.ModifiedFEMBaseFunction;
import ua.lviv.franko.mlsFunctions.oneDimension.Particle;
import ua.lviv.franko.mlsFunctions.oneDimension.weightFunctions.CubicSpline;
import ua.lviv.franko.nonequidistant.FEMShapeFunctions1D.FEMShapeFunction1D1Pow;
import ua.lviv.franko.visualisation.NewOkCancelDialog;

public class SolverCoupled {

	protected ArrayList<ICalculate>	baseFunctionsList;
	protected double				left;
	protected double				right;
	protected EquationParameters	eParams;

	protected HashMatrixLUStyle			mainMatrix;
	protected double[]				b;
	public static int				numberOfQuadraturePoints	= 20;

	public SolverCoupled(ArrayList<ICalculate> baseFunctionsList, double left, double right, EquationParameters eParams) {
		this.eParams = eParams;
		this.left = left;
		this.right = right;
		this.baseFunctionsList = baseFunctionsList;
	}

	public void generate() {
		Quadrature gau = new GaussQuadrature(numberOfQuadraturePoints);
		b = new double[baseFunctionsList.size()];
		mainMatrix = new HashMatrixLUStyle(RowType.ARRAY_TYPE, baseFunctionsList.size());
		for (int i = 0; i < b.length; i++) {
			final ICalculate phi1 = baseFunctionsList.get(i);
			double[] bounds = phi1.getBounds();
			bounds[0] = Math.max(left, bounds[0]);
			bounds[1] = Math.min(right, bounds[1]);
			for (int j = 0; j < b.length; j++) {
				final ICalculate phi2 = baseFunctionsList.get(j);
				// bounds[0] = Math.max(left, Math.max(phi1.getBounds()[0],
				// phi2.getBounds()[0]));
				// bounds[1] = Math.min(right, Math.min(phi1.getBounds()[1],
				// phi2.getBounds()[1]));
				double value = gau.integrate(bounds[0], bounds[1], new IFunction1D() {

					@Override
					public double calculate(double x) {
						double[] r1 = phi2.calculateAll(x);
						double[] r2 = phi1.calculateAll(x);
						return eParams.p.calculate(x) * r1[1] * r2[1] + eParams.w.calculate(x) * r1[1] * r2[0] + eParams.q.calculate(x) * r1[0] * r2[0];
					}

					@Override
					public double calculateDerivaty(double x) {

						return 0;
					}
				});
				mainMatrix.setValue(i, j, value);
			}
			// bounds = phi1.getBounds();
			b[i] = gau.integrate(bounds[0], bounds[1], new IFunction1D() {

				@Override
				public double calculate(double x) {
					double r1 = phi1.calculate(x);
					return eParams.f.calculate(x) * r1;
				}

				@Override
				public double calculateDerivaty(double x) {
					return 0.0;
				}
			});
		}
	}

	public double[] solve() {

		mainMatrix.getJRowValues(0).clear();
		mainMatrix.getJRowIndex(0).clear();

		mainMatrix.setValue(0, 0, 1);
		b[0] = 0.0;

		mainMatrix.getJRowValues(mainMatrix.getN() - 1).clear();
		mainMatrix.getJRowIndex(mainMatrix.getN() - 1).clear();
		mainMatrix.setValue(mainMatrix.getN() - 1, mainMatrix.getN() - 1, 1);
		b[b.length - 1] = 0.0;

		LUFactorization lu = new LUFactorization(this.mainMatrix);
		return lu.calculate(b);
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		ArrayList<IFEM1D> fem = new ArrayList<IFEM1D>();
		// fem.add(new FEMShapeFunction1D1(0.0, 0.005));
		// fem.add(new FEMShapeFunction1D1(1.0, 0.005));

		final ArrayList<ICalculate> baseFuncions = new ArrayList<ICalculate>();
		// fem.add(new FEMShapeFunction1D1(0.0, 0.1));
		// fem.add(new FEMShapeFunction1D1(0.1, 0.1));
		// fem.add(new FEMShapeFunction1D1(0.2, 0.1));
		// fem.add(new FEMShapeFunction1D1(0.3, 0.1));
		// fem.add(new FEMShapeFunction1D1(0.4, 0.1));
		int fn = 4;
		double divPoint = 0.85;
		double fh = (divPoint) / (double) fn;
		if (fn == 0) {
			fh = 0.01;
			divPoint = 0.0;
		}
		for (int i = 0; i < fn; i++) {
			// fem.add(new FEMShapeFunction1D1(fh*i, fh));
			fem.add(new FEMShapeFunction1D1Pow(fh * i, fh, fh));
		}

		double a = divPoint, b = 1;
		int nn = 20;
		ArrayList<Particle> ls1D = new ArrayList<Particle>();
		int PolinomDegree = 1;
		double h = (b - a) / (double) nn;
		double rho = 2.33 * h;
		for (int i = 0; i <= nn; i++) {
			ls1D.add(new Particle(a + h * i, rho));
		}

		a = 0.0;

		for (Iterator<IFEM1D> iterator = fem.iterator(); iterator.hasNext();) {
			ICalculate iCalculate = (ICalculate) iterator.next();
			baseFuncions.add(iCalculate);
		}
		fem.add(new FEMShapeFunction1D1Pow(divPoint, fh, h));
		// fem.add(new FEMShapeFunction1D1Pow(h, h,h));
		// fem.add(new FEMShapeFunction1D1Pow(2*h, h,h));
		// fem.add(new FEMShapeFunction1D1Pow(3*h, h,h));
		// fem.add(new FEMShapeFunction1D1Pow(4*h, h,h));
		// fem.add(new FEMShapeFunction1D1Pow(5*h, h,h));
		fem.add(new FEMShapeFunction1D1Pow(1.0, h, fh));

		CubicSpline cb = new CubicSpline();

		// ModifiedFEMBaseFunctiion mf = new
		// ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(0), cb, fem,
		// fem.get(fem.size()-2));
		ModifiedFEMBaseFunction mf = new ModifiedFEMBaseFunction(PolinomDegree, ls1D.get(0), cb, fem, fem.get(0));
		mf.setNodeList(ls1D);
		baseFuncions.add(mf);
		// mf = new ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(1), cb,
		// fem, fem.get(1));
		// mf.setNodeList(ls1D);
		// baseFuncions.add(mf);
		// mf = new ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(2), cb,
		// fem, fem.get(2));
		// mf.setNodeList(ls1D);
		// baseFuncions.add(mf);
		// mf = new ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(3), cb,
		// fem, fem.get(3));
		// mf.setNodeList(ls1D);
		// baseFuncions.add(mf);
		// mf = new ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(4), cb,
		// fem, fem.get(4));
		// mf.setNodeList(ls1D);
		// baseFuncions.add(mf);
		// mf = new ModifiedFEMBaseFunctiion(PolinomDegree, ls1D.get(5), cb,
		// fem, fem.get(5));
		// mf.setNodeList(ls1D);
		// baseFuncions.add(mf);
		for (int i = 1; i < ls1D.size() - 1; i++) {
			Particle particle = ls1D.get(i);
			BaseFunctionCoupled jmf1 = new BaseFunctionCoupled(PolinomDegree, particle, cb, fem);
			jmf1.setNodeList(ls1D);
			baseFuncions.add(jmf1);
		}
		mf = new ModifiedFEMBaseFunction(PolinomDegree, ls1D.get(ls1D.size() - 1), cb, fem, fem.get(fem.size() - 1));
		mf.setNodeList(ls1D);
		baseFuncions.add(mf);

		XYSeriesCollection col = new XYSeriesCollection();
		// for (Iterator<ICalculate> iterator = baseFuncions.iterator();
		// iterator.hasNext();) {
		//
		//
		// XYSeries s = new XYSeries("1", true, true);
		// ICalculate jmf1 = iterator.next();
		// System.out.println("x = [");
		// double xxx= a;
		// while(xxx<=b+0.0001){
		// double tmp1 = jmf1.calculateAll(xxx)[0];
		// s.add(xxx,tmp1);
		// System.out.println(xxx+"\t"+tmp1);
		// xxx+=0.002;
		// }
		// System.out.println("];\nplot(x(:,1),x(:,2));");
		// col.addSeries(s);
		// }

		// XYSeries s = new XYSeries("1", true, true);
		// ICalculate jmf1 = baseFuncions.get(0);
		// double xxx= a;
		// while(xxx<=b+0.0001){
		// double tmp1 = jmf1.calculateAll(xxx)[0];
		// s.add(xxx,tmp1);
		// xxx+=0.002;
		// }
		// col.addSeries(s);
		//
		// s = new XYSeries("1", true, true);
		// jmf1 = baseFuncions.get(baseFuncions.size()-1);
		// xxx= a;
		// while(xxx<=b+0.0001){
		// double tmp1 = jmf1.calculateAll(xxx)[0];
		// s.add(xxx,tmp1);
		// xxx+=0.002;
		// }
		// col.addSeries(s);

		NewOkCancelDialog d = new NewOkCancelDialog(null, false, col);
		// d.setLocationRelativeTo(null);
		// d.setVisible(true);

		EquationParameters eParams = new EquationParameters(new ConstantFunction1D(1.0), new ConstantFunction1D(10.0), new ConstantFunction1D(10.0),
				new FunctionWithoutDerivatives1D() {

					@Override
					public double calculate(double x) {
						return 25.0 * Math.sin(x / 4.0);
					}
				});

		SolverCoupled solver = new SolverCoupled(baseFuncions, a, b, eParams);
		solver.generate();

		final double[] res = solver.solve();

		XYSeries s2 = new XYSeries("result", true, true);
		System.out.println("x = [");
		for (double x = 0; x < b + 0.001; x += 0.001) {
			double val = 0.0;
			for (int k = 0; k < res.length; k++) {
				ICalculate func = baseFuncions.get(k);
				val += func.calculate(x) * res[k];
			}
			System.out.println(x + "\t" + val);
			s2.add(x, val);
		}
		System.out.println("];\nplot(x(:,1),x(:,2));");
		col = new XYSeriesCollection();
		col.addSeries(s2);
		d = new NewOkCancelDialog(null, false, col);
		d.setLocationRelativeTo(null);
		d.setVisible(true);

		Quadrature gau = new GaussQuadrature(100);
		System.out.print(gau.integrate(a, b, new IFunction1D() {

			@Override
			public double calculate(double x) {
				double tmp = 0.0;
				for (int j = 0; j < res.length; j++) {
					tmp += baseFuncions.get(j).calculate(x) * res[j];
				}

				double nu = 0.01;
				double A = 1;
				double aux = Math.PI * (A * A + nu * nu * Math.PI * Math.PI);
				double e = Math.exp(A / nu);
				double c1 = (-aux + A * (e + 1)) / (aux * (e - 1));
				double c2 = (aux - 2 * A) / (aux * (e - 1));
				double r = c1 + c2 * Math.exp(A * x / nu) + nu * Math.PI * (Math.sin(Math.PI * x) - A * Math.cos(Math.PI * x) / (nu * Math.PI)) / aux;

				tmp -= (r);
				return tmp * tmp;
			}

			@Override
			public double calculateDerivaty(double x) {

				return 0;
			}
		}));

	}

}
